Load the student scores for the test - here we load the ETH Zurich test data, downloaded from https://pontifex.ethz.ch/s21t5/rate/
head(test_scores_pre)
## # A tibble: 6 x 38
## year class A1 A2 A3 A4 A5 A6 A7 A8 A9 A10 A11
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2017 s21t-~ 1 0 1 1 1 0 1 0 1 0 1
## 2 2017 s21t-~ 1 0 1 1 1 1 0 1 1 2 1
## 3 2017 s21t-~ 1 0 0 0 1 1 1 0 1 1 1
## 4 2017 s21t-~ 1 0 1 1 1 1 1 1 1 1 1
## 5 2017 s21t-~ 1 0 1 0 2 0 1 0 1 0 2
## 6 2017 s21t-~ 0 1 0 0 1 2 0 2 2 2 1
## # ... with 25 more variables: A12 <dbl>, A13 <dbl>, A14 <dbl>, A15 <dbl>,
## # A16 <dbl>, A17 <dbl>, A18 <dbl>, A19 <dbl>, A20 <dbl>, A21 <dbl>,
## # A22 <dbl>, A23 <dbl>, A24 <dbl>, A25 <dbl>, A26 <dbl>, A27 <dbl>,
## # A28 <dbl>, A29 <dbl>, A30 <dbl>, A31 <dbl>, A32 <dbl>, A33 <dbl>,
## # A34 <dbl>, A35 <dbl>, A36 <dbl>
head(test_scores_post)
## # A tibble: 6 x 34
## year class B1 B2 B3 B4 B5 B6 B7 B8 B9 B10 B11
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2019 s21t-~ 1 1 1 1 1 1 1 1 1 1 0
## 2 2019 s21t-~ 2 1 0 0 0 1 2 1 2 0 1
## 3 2019 s21t-~ 1 1 0 1 1 1 1 1 1 1 0
## 4 2019 s21t-~ 1 1 1 0 1 0 0 0 1 1 1
## 5 2019 s21t-~ 1 0 1 0 0 1 1 0 0 1 0
## 6 2019 s21t-~ 1 0 0 0 1 1 0 0 1 1 0
## # ... with 21 more variables: B12 <dbl>, B13 <dbl>, B14 <dbl>, B15 <dbl>,
## # B16 <dbl>, B17 <dbl>, B18 <dbl>, B19 <dbl>, B20 <dbl>, B21 <dbl>,
## # B22 <dbl>, B23 <dbl>, B24 <dbl>, B25 <dbl>, B26 <dbl>, B27 <dbl>,
## # B28 <dbl>, B29 <dbl>, B30 <dbl>, B31 <dbl>, B32 <dbl>
The number of responses from each cohort:
test_scores %>%
group_by(year) %>%
tally() %>%
gt() %>%
data_color(
columns = c("n"),
colors = scales::col_numeric(palette = c("Blues"), domain = NULL)
)
## Warning: The `.dots` argument of `group_by()` is deprecated as of dplyr 1.0.0.
| year | n |
|---|---|
| 2017 | 1678 |
| 2018 | 1746 |
| 2019 | 1995 |
| 2020 | 2188 |
| 2021 | 2041 |
Mean and standard deviation for each item:
test_scores %>%
select(-class) %>%
group_by(year) %>%
skim_without_charts() %>%
select(-contains("character."), -contains("numeric.p"), -skim_type) %>%
group_by(year) %>%
gt() %>%
fmt_number(columns = contains("numeric"), decimals = 3) %>%
data_color(
columns = c("numeric.mean"),
colors = scales::col_numeric(palette = c("Greens"), domain = NULL)
) %>%
cols_label(
numeric.mean = "Mean",
numeric.sd = "SD"
)
| skim_variable | n_missing | complete_rate | Mean | SD |
|---|---|---|---|---|
| 2017 | ||||
| test_version | 0 | 1.0000000 | NA | NA |
| A1_B1 | 15 | 0.9910608 | 0.960 | 0.195 |
| A2_B0 | 131 | 0.9219309 | 0.333 | 0.471 |
| A3_B2 | 77 | 0.9541120 | 0.655 | 0.476 |
| A4_B3 | 27 | 0.9839094 | 0.655 | 0.476 |
| A5_B4 | 544 | 0.6758045 | 0.692 | 0.462 |
| A6_B5 | 362 | 0.7842670 | 0.901 | 0.298 |
| A7_B6 | 19 | 0.9886770 | 0.716 | 0.451 |
| A8_B7 | 311 | 0.8146603 | 0.583 | 0.493 |
| A9_B8 | 613 | 0.6346841 | 0.731 | 0.444 |
| A10_B9 | 452 | 0.7306317 | 0.745 | 0.436 |
| A11_B10 | 352 | 0.7902265 | 0.729 | 0.445 |
| A12_B0 | 216 | 0.8712753 | 0.784 | 0.412 |
| A13_B0 | 211 | 0.8742551 | 0.429 | 0.495 |
| A14_B12 | 237 | 0.8587604 | 0.653 | 0.476 |
| A15_B13 | 269 | 0.8396901 | 0.551 | 0.498 |
| A16_B14 | 427 | 0.7455304 | 0.249 | 0.433 |
| A17_B15 | 154 | 0.9082241 | 0.649 | 0.477 |
| A18_B16 | 468 | 0.7210965 | 0.529 | 0.499 |
| A19_B0 | 566 | 0.6626937 | 0.533 | 0.499 |
| A20_B17 | 108 | 0.9356377 | 0.776 | 0.417 |
| A21_B18 | 408 | 0.7568534 | 0.674 | 0.469 |
| A22_B19 | 477 | 0.7157330 | 0.723 | 0.448 |
| A23_B20 | 492 | 0.7067938 | 0.584 | 0.493 |
| A24_B21 | 326 | 0.8057211 | 0.730 | 0.444 |
| A25_B0 | 78 | 0.9535161 | 0.568 | 0.495 |
| A26_B22 | 250 | 0.8510131 | 0.938 | 0.241 |
| A27_B23 | 590 | 0.6483909 | 0.603 | 0.490 |
| A28_B24 | 647 | 0.6144219 | 0.686 | 0.464 |
| A29_B25 | 373 | 0.7777116 | 0.624 | 0.485 |
| A30_B26 | 145 | 0.9135876 | 0.830 | 0.376 |
| A31_B27 | 214 | 0.8724672 | 0.432 | 0.495 |
| A32_B28 | 466 | 0.7222884 | 0.315 | 0.465 |
| A33_B29 | 57 | 0.9660310 | 0.808 | 0.394 |
| A34_B30 | 497 | 0.7038141 | 0.826 | 0.379 |
| A35_B0 | 886 | 0.4719905 | 0.688 | 0.464 |
| A36_B0 | 830 | 0.5053635 | 0.436 | 0.496 |
| A0_B11 | 1678 | 0.0000000 | NaN | NA |
| A0_B31 | 1678 | 0.0000000 | NaN | NA |
| A0_B32 | 1678 | 0.0000000 | NaN | NA |
| 2018 | ||||
| test_version | 0 | 1.0000000 | NA | NA |
| A1_B1 | 16 | 0.9908362 | 0.964 | 0.186 |
| A2_B0 | 97 | 0.9444444 | 0.300 | 0.458 |
| A3_B2 | 37 | 0.9788087 | 0.658 | 0.475 |
| A4_B3 | 23 | 0.9868270 | 0.634 | 0.482 |
| A5_B4 | 441 | 0.7474227 | 0.659 | 0.474 |
| A6_B5 | 269 | 0.8459336 | 0.870 | 0.336 |
| A7_B6 | 13 | 0.9925544 | 0.711 | 0.453 |
| A8_B7 | 199 | 0.8860252 | 0.561 | 0.496 |
| A9_B8 | 549 | 0.6855670 | 0.673 | 0.469 |
| A10_B9 | 367 | 0.7898053 | 0.687 | 0.464 |
| A11_B10 | 293 | 0.8321879 | 0.700 | 0.458 |
| A12_B0 | 127 | 0.9272623 | 0.770 | 0.421 |
| A13_B0 | 188 | 0.8923253 | 0.401 | 0.490 |
| A14_B12 | 192 | 0.8900344 | 0.645 | 0.479 |
| A15_B13 | 220 | 0.8739977 | 0.530 | 0.499 |
| A16_B14 | 341 | 0.8046964 | 0.250 | 0.433 |
| A17_B15 | 120 | 0.9312715 | 0.648 | 0.478 |
| A18_B16 | 372 | 0.7869416 | 0.463 | 0.499 |
| A19_B0 | 468 | 0.7319588 | 0.468 | 0.499 |
| A20_B17 | 85 | 0.9513173 | 0.762 | 0.426 |
| A21_B18 | 336 | 0.8075601 | 0.646 | 0.478 |
| A22_B19 | 420 | 0.7594502 | 0.682 | 0.466 |
| A23_B20 | 397 | 0.7726231 | 0.558 | 0.497 |
| A24_B21 | 260 | 0.8510882 | 0.678 | 0.467 |
| A25_B0 | 46 | 0.9736541 | 0.753 | 0.431 |
| A26_B22 | 230 | 0.8682703 | 0.929 | 0.256 |
| A27_B23 | 531 | 0.6958763 | 0.576 | 0.494 |
| A28_B24 | 529 | 0.6970218 | 0.655 | 0.476 |
| A29_B25 | 319 | 0.8172967 | 0.633 | 0.482 |
| A30_B26 | 121 | 0.9306987 | 0.786 | 0.410 |
| A31_B27 | 166 | 0.9049255 | 0.411 | 0.492 |
| A32_B28 | 436 | 0.7502864 | 0.272 | 0.445 |
| A33_B29 | 63 | 0.9639175 | 0.791 | 0.407 |
| A34_B30 | 425 | 0.7565865 | 0.815 | 0.389 |
| A35_B0 | 797 | 0.5435281 | 0.612 | 0.488 |
| A36_B0 | 741 | 0.5756014 | 0.397 | 0.490 |
| A0_B11 | 1746 | 0.0000000 | NaN | NA |
| A0_B31 | 1746 | 0.0000000 | NaN | NA |
| A0_B32 | 1746 | 0.0000000 | NaN | NA |
| 2019 | ||||
| test_version | 0 | 1.0000000 | NA | NA |
| A1_B1 | 12 | 0.9939850 | 0.963 | 0.188 |
| A2_B0 | 1995 | 0.0000000 | NaN | NA |
| A3_B2 | 29 | 0.9854637 | 0.678 | 0.467 |
| A4_B3 | 35 | 0.9824561 | 0.655 | 0.476 |
| A5_B4 | 513 | 0.7428571 | 0.661 | 0.473 |
| A6_B5 | 296 | 0.8516291 | 0.875 | 0.331 |
| A7_B6 | 25 | 0.9874687 | 0.701 | 0.458 |
| A8_B7 | 195 | 0.9022556 | 0.582 | 0.493 |
| A9_B8 | 628 | 0.6852130 | 0.670 | 0.470 |
| A10_B9 | 443 | 0.7779449 | 0.713 | 0.453 |
| A11_B10 | 299 | 0.8501253 | 0.685 | 0.465 |
| A12_B0 | 1995 | 0.0000000 | NaN | NA |
| A13_B0 | 1995 | 0.0000000 | NaN | NA |
| A14_B12 | 229 | 0.8852130 | 0.664 | 0.473 |
| A15_B13 | 278 | 0.8606516 | 0.550 | 0.498 |
| A16_B14 | 372 | 0.8135338 | 0.251 | 0.434 |
| A17_B15 | 133 | 0.9333333 | 0.641 | 0.480 |
| A18_B16 | 487 | 0.7558897 | 0.366 | 0.482 |
| A19_B0 | 1995 | 0.0000000 | NaN | NA |
| A20_B17 | 101 | 0.9493734 | 0.778 | 0.416 |
| A21_B18 | 368 | 0.8155388 | 0.661 | 0.473 |
| A22_B19 | 453 | 0.7729323 | 0.669 | 0.471 |
| A23_B20 | 439 | 0.7799499 | 0.540 | 0.499 |
| A24_B21 | 268 | 0.8656642 | 0.690 | 0.463 |
| A25_B0 | 1995 | 0.0000000 | NaN | NA |
| A26_B22 | 261 | 0.8691729 | 0.927 | 0.261 |
| A27_B23 | 580 | 0.7092732 | 0.595 | 0.491 |
| A28_B24 | 600 | 0.6992481 | 0.649 | 0.478 |
| A29_B25 | 344 | 0.8275689 | 0.614 | 0.487 |
| A30_B26 | 240 | 0.8796992 | 0.612 | 0.487 |
| A31_B27 | 192 | 0.9037594 | 0.438 | 0.496 |
| A32_B28 | 446 | 0.7764411 | 0.260 | 0.439 |
| A33_B29 | 69 | 0.9654135 | 0.799 | 0.401 |
| A34_B30 | 420 | 0.7894737 | 0.770 | 0.421 |
| A35_B0 | 1995 | 0.0000000 | NaN | NA |
| A36_B0 | 1995 | 0.0000000 | NaN | NA |
| A0_B11 | 124 | 0.9378446 | 0.607 | 0.489 |
| A0_B31 | 618 | 0.6902256 | 0.502 | 0.500 |
| A0_B32 | 1179 | 0.4090226 | 0.488 | 0.500 |
| 2020 | ||||
| test_version | 0 | 1.0000000 | NA | NA |
| A1_B1 | 13 | 0.9940585 | 0.963 | 0.188 |
| A2_B0 | 2188 | 0.0000000 | NaN | NA |
| A3_B2 | 48 | 0.9780622 | 0.671 | 0.470 |
| A4_B3 | 31 | 0.9858318 | 0.632 | 0.482 |
| A5_B4 | 592 | 0.7294333 | 0.670 | 0.470 |
| A6_B5 | 299 | 0.8633455 | 0.886 | 0.318 |
| A7_B6 | 20 | 0.9908592 | 0.711 | 0.453 |
| A8_B7 | 283 | 0.8706581 | 0.691 | 0.462 |
| A9_B8 | 661 | 0.6978976 | 0.623 | 0.485 |
| A10_B9 | 481 | 0.7801645 | 0.705 | 0.456 |
| A11_B10 | 354 | 0.8382084 | 0.710 | 0.454 |
| A12_B0 | 2188 | 0.0000000 | NaN | NA |
| A13_B0 | 2188 | 0.0000000 | NaN | NA |
| A14_B12 | 253 | 0.8843693 | 0.613 | 0.487 |
| A15_B13 | 279 | 0.8724863 | 0.524 | 0.500 |
| A16_B14 | 410 | 0.8126143 | 0.243 | 0.429 |
| A17_B15 | 139 | 0.9364717 | 0.632 | 0.483 |
| A18_B16 | 492 | 0.7751371 | 0.351 | 0.478 |
| A19_B0 | 2188 | 0.0000000 | NaN | NA |
| A20_B17 | 100 | 0.9542962 | 0.741 | 0.438 |
| A21_B18 | 401 | 0.8167276 | 0.632 | 0.482 |
| A22_B19 | 481 | 0.7801645 | 0.664 | 0.473 |
| A23_B20 | 512 | 0.7659963 | 0.546 | 0.498 |
| A24_B21 | 313 | 0.8569470 | 0.659 | 0.474 |
| A25_B0 | 2188 | 0.0000000 | NaN | NA |
| A26_B22 | 276 | 0.8738574 | 0.929 | 0.256 |
| A27_B23 | 635 | 0.7097806 | 0.558 | 0.497 |
| A28_B24 | 710 | 0.6755027 | 0.616 | 0.487 |
| A29_B25 | 409 | 0.8130713 | 0.604 | 0.489 |
| A30_B26 | 288 | 0.8683729 | 0.615 | 0.487 |
| A31_B27 | 268 | 0.8775137 | 0.461 | 0.499 |
| A32_B28 | 529 | 0.7582267 | 0.274 | 0.446 |
| A33_B29 | 76 | 0.9652651 | 0.763 | 0.425 |
| A34_B30 | 464 | 0.7879342 | 0.770 | 0.421 |
| A35_B0 | 2188 | 0.0000000 | NaN | NA |
| A36_B0 | 2188 | 0.0000000 | NaN | NA |
| A0_B11 | 123 | 0.9437843 | 0.630 | 0.483 |
| A0_B31 | 676 | 0.6910420 | 0.472 | 0.499 |
| A0_B32 | 1266 | 0.4213894 | 0.507 | 0.500 |
| 2021 | ||||
| test_version | 0 | 1.0000000 | NA | NA |
| A1_B1 | 16 | 0.9921607 | 0.963 | 0.189 |
| A2_B0 | 2041 | 0.0000000 | NaN | NA |
| A3_B2 | 35 | 0.9828515 | 0.680 | 0.466 |
| A4_B3 | 42 | 0.9794219 | 0.617 | 0.486 |
| A5_B4 | 513 | 0.7486526 | 0.657 | 0.475 |
| A6_B5 | 297 | 0.8544831 | 0.872 | 0.334 |
| A7_B6 | 16 | 0.9921607 | 0.687 | 0.464 |
| A8_B7 | 261 | 0.8721215 | 0.669 | 0.471 |
| A9_B8 | 606 | 0.7030867 | 0.606 | 0.489 |
| A10_B9 | 474 | 0.7677609 | 0.690 | 0.462 |
| A11_B10 | 360 | 0.8236159 | 0.694 | 0.461 |
| A12_B0 | 2041 | 0.0000000 | NaN | NA |
| A13_B0 | 2041 | 0.0000000 | NaN | NA |
| A14_B12 | 260 | 0.8726115 | 0.616 | 0.487 |
| A15_B13 | 267 | 0.8691818 | 0.525 | 0.500 |
| A16_B14 | 404 | 0.8020578 | 0.246 | 0.431 |
| A17_B15 | 122 | 0.9402254 | 0.611 | 0.488 |
| A18_B16 | 449 | 0.7800098 | 0.358 | 0.480 |
| A19_B0 | 2041 | 0.0000000 | NaN | NA |
| A20_B17 | 105 | 0.9485546 | 0.750 | 0.433 |
| A21_B18 | 405 | 0.8015679 | 0.639 | 0.481 |
| A22_B19 | 479 | 0.7653111 | 0.653 | 0.476 |
| A23_B20 | 461 | 0.7741303 | 0.549 | 0.498 |
| A24_B21 | 325 | 0.8407643 | 0.681 | 0.466 |
| A25_B0 | 2041 | 0.0000000 | NaN | NA |
| A26_B22 | 259 | 0.8731014 | 0.919 | 0.273 |
| A27_B23 | 609 | 0.7016169 | 0.585 | 0.493 |
| A28_B24 | 646 | 0.6834885 | 0.602 | 0.490 |
| A29_B25 | 387 | 0.8103871 | 0.619 | 0.486 |
| A30_B26 | 290 | 0.8579128 | 0.613 | 0.487 |
| A31_B27 | 280 | 0.8628123 | 0.469 | 0.499 |
| A32_B28 | 485 | 0.7623714 | 0.269 | 0.443 |
| A33_B29 | 79 | 0.9612935 | 0.783 | 0.412 |
| A34_B30 | 451 | 0.7790299 | 0.771 | 0.420 |
| A35_B0 | 2041 | 0.0000000 | NaN | NA |
| A36_B0 | 2041 | 0.0000000 | NaN | NA |
| A0_B11 | 106 | 0.9480647 | 0.613 | 0.487 |
| A0_B31 | 519 | 0.7457129 | 0.604 | 0.489 |
| A0_B32 | 1154 | 0.4345909 | 0.488 | 0.500 |
Before applying IRT, we should check that the data satisfies the assumptions needed by the model. In particular, to use a 1-dimensional IRT model, we should have some evidence of unidimensionality in the test scores.
This plot shows the correlations between scores on each pair of items – note that it is restricted to only those items that appear on both versions of the test, since the plotting package did not deal well with missing data:
item_scores <- test_scores %>%
select(-class, -year, -test_version)
item_scores_unchanged_only <- item_scores %>%
select(!contains("B0")) %>% select(!contains("A0"))
cor_ci <- psych::corCi(item_scores_unchanged_only, plot = FALSE)
psych::cor.plot.upperLowerCi(cor_ci)
There are a few correlations that are not significantly different from 0:
cor_ci$ci %>%
as_tibble(rownames = "corr") %>%
filter(p > 0.05) %>%
arrange(-p) %>%
select(-contains(".e")) %>%
gt() %>%
fmt_number(columns = 2:4, decimals = 3)
| corr | lower | upper | p |
|---|---|---|---|
| A1_B1-A29_B | −0.016 | 0.032 | 0.535 |
| A29_B-A31_B | −0.040 | 0.013 | 0.305 |
| A24_B-A29_B | −0.008 | 0.039 | 0.203 |
Here we redo the correlation calculations with all the items, and check that there are still few cases where the correlations close to 0:
cor_ci2 <- psych::corCi(item_scores, plot = FALSE)
cor_ci2$ci %>%
as_tibble(rownames = "corr") %>%
filter(p > 0.05) %>%
arrange(-p) %>%
select(-contains(".e")) %>%
gt() %>%
fmt_number(columns = 2:4, decimals = 3)
| corr | lower | upper | p |
|---|---|---|---|
| A1_B1-A29_B | −0.014 | 0.029 | 0.494 |
| A2_B0-A34_B | −0.019 | 0.063 | 0.284 |
| A24_B-A29_B | −0.009 | 0.040 | 0.215 |
| A29_B-A31_B | −0.037 | 0.008 | 0.206 |
| A1_B1-A2_B0 | −0.011 | 0.061 | 0.176 |
| A2_B0-A25_B | −0.010 | 0.070 | 0.142 |
| A1_B1-A12_B | −0.006 | 0.076 | 0.096 |
| A25_B-A35_B | −0.004 | 0.098 | 0.069 |
| A12_B-A34_B | −0.002 | 0.093 | 0.059 |
The overall picture is that the item scores are well correlated with each other.
Here we again focus on the subset of items that appeared in both tests.
structure <- check_factorstructure(item_scores_unchanged_only)
n <- n_factors(item_scores_unchanged_only)
The choice of 1 dimensions is supported by 6 (26.09%) methods out of 23 (Acceleration factor, R2, VSS complexity 1, Velicer’s MAP, TLI, RMSEA).
plot(n)
summary(n) %>% gt()
| n_Factors | n_Methods |
|---|---|
| 1 | 6 |
| 2 | 2 |
| 3 | 1 |
| 4 | 6 |
| 7 | 1 |
| 11 | 2 |
| 19 | 1 |
| 21 | 1 |
| 28 | 3 |
#n %>% tibble() %>% gt()
fa.parallel(item_scores_unchanged_only, fa = "fa")
## Parallel analysis suggests that the number of factors = 8 and the number of components = NA
item_scores_unchanged_only_and_no_na <- item_scores_unchanged_only %>%
mutate(
across(everything(), ~replace_na(.x, 0))
)
fitfact <- factanal(item_scores_unchanged_only_and_no_na, factors = 1, rotation = "varimax")
print(fitfact, digits = 2, cutoff = 0.3, sort = TRUE)
##
## Call:
## factanal(x = item_scores_unchanged_only_and_no_na, factors = 1, rotation = "varimax")
##
## Uniquenesses:
## A1_B1 A3_B2 A4_B3 A5_B4 A6_B5 A7_B6 A8_B7 A9_B8 A10_B9 A11_B10
## 0.96 0.77 0.80 0.73 0.77 0.77 0.79 0.62 0.66 0.69
## A14_B12 A15_B13 A16_B14 A17_B15 A18_B16 A20_B17 A21_B18 A22_B19 A23_B20 A24_B21
## 0.72 0.88 0.82 0.78 0.84 0.65 0.59 0.61 0.62 0.74
## A26_B22 A27_B23 A28_B24 A29_B25 A30_B26 A31_B27 A32_B28 A33_B29 A34_B30
## 0.70 0.64 0.64 0.87 0.82 0.77 0.80 0.88 0.88
##
## Loadings:
## [1] 0.52 0.62 0.59 0.56 0.53 0.59 0.64 0.62 0.61 0.51 0.55 0.60 0.60 0.48
## [16] 0.45 0.48 0.48 0.46 0.35 0.42 0.46 0.40 0.36 0.42 0.48 0.45 0.35 0.34
##
## Factor1
## SS loadings 7.19
## Proportion Var 0.25
##
## Test of the hypothesis that 1 factor is sufficient.
## The chi square statistic is 5073.36 on 377 degrees of freedom.
## The p-value is 0
load <- tidy(fitfact)
ggplot(load, aes(x = fl1, y = 0)) +
geom_point() +
geom_label_repel(aes(label = paste0("A", rownames(load))), show.legend = FALSE) +
labs(x = "Factor 1", y = NULL,
title = "Standardised Loadings",
subtitle = "Based upon correlation matrix") +
theme_minimal()
## Warning: ggrepel: 10 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
# To proceed with the IRT analysis, comment out the following line before knitting
#knitr::knit_exit()
We can fit a Multidimensional Item Response Theory (mirt) model. From the function definition:
mirt fits a maximum likelihood (or maximum a posteriori) factor analysis model to any mixture of dichotomous and polytomous data under the item response theory paradigm using either Cai's (2010) Metropolis-Hastings Robbins-Monro (MHRM) algorithm.
The process is to first fit the model, and save the result as a model object that we can then parse to get tabular or visual displays of the model that we might want. When fitting the model, we have the option to specify a few arguments, which then get interpreted as parameters to be passed to the model.
fit_2pl <- mirt(
data = item_scores, # just the columns with question scores
#removeEmptyRows = TRUE,
model = 1, # number of factors to extract
itemtype = "2PL", # 2 parameter logistic model
SE = TRUE # estimate standard errors
)
##
Iteration: 1, Log-Lik: -153078.292, Max-Change: 3.97531
Iteration: 2, Log-Lik: -144683.865, Max-Change: 3.18540
Iteration: 3, Log-Lik: -143129.397, Max-Change: 0.51288
Iteration: 4, Log-Lik: -142627.917, Max-Change: 0.56576
Iteration: 5, Log-Lik: -142393.297, Max-Change: 0.16584
Iteration: 6, Log-Lik: -142217.904, Max-Change: 0.23612
Iteration: 7, Log-Lik: -142107.527, Max-Change: 0.10023
Iteration: 8, Log-Lik: -142014.839, Max-Change: 0.13562
Iteration: 9, Log-Lik: -141948.534, Max-Change: 0.06224
Iteration: 10, Log-Lik: -141891.498, Max-Change: 0.09381
Iteration: 11, Log-Lik: -141852.472, Max-Change: 0.06481
Iteration: 12, Log-Lik: -141816.748, Max-Change: 0.06046
Iteration: 13, Log-Lik: -141792.864, Max-Change: 0.03712
Iteration: 14, Log-Lik: -141771.799, Max-Change: 0.04925
Iteration: 15, Log-Lik: -141757.902, Max-Change: 0.03397
Iteration: 16, Log-Lik: -141746.237, Max-Change: 0.02459
Iteration: 17, Log-Lik: -141737.341, Max-Change: 0.02003
Iteration: 18, Log-Lik: -141730.332, Max-Change: 0.01591
Iteration: 19, Log-Lik: -141724.992, Max-Change: 0.01455
Iteration: 20, Log-Lik: -141720.793, Max-Change: 0.01250
Iteration: 21, Log-Lik: -141717.574, Max-Change: 0.01163
Iteration: 22, Log-Lik: -141711.914, Max-Change: 0.01239
Iteration: 23, Log-Lik: -141710.375, Max-Change: 0.00837
Iteration: 24, Log-Lik: -141709.314, Max-Change: 0.00617
Iteration: 25, Log-Lik: -141707.449, Max-Change: 0.00420
Iteration: 26, Log-Lik: -141707.143, Max-Change: 0.00409
Iteration: 27, Log-Lik: -141706.907, Max-Change: 0.00287
Iteration: 28, Log-Lik: -141706.461, Max-Change: 0.00266
Iteration: 29, Log-Lik: -141706.391, Max-Change: 0.00172
Iteration: 30, Log-Lik: -141706.338, Max-Change: 0.00152
Iteration: 31, Log-Lik: -141706.204, Max-Change: 0.00016
Iteration: 32, Log-Lik: -141706.203, Max-Change: 0.00012
Iteration: 33, Log-Lik: -141706.202, Max-Change: 0.00012
Iteration: 34, Log-Lik: -141706.198, Max-Change: 0.00045
Iteration: 35, Log-Lik: -141706.197, Max-Change: 0.00009
##
## Calculating information matrix...
We then compute factor score estimates and augment the existing data frame with these estimates, to keep everything in one place. To do the estimation, we use the fscores() function from the mirt package which takes in a computed model object and computes factor score estimates according to the method specified. We will use the EAP method for factor score estimation, which is the “expected a-posteriori” method, the default. We specify it explicitly below, but the results would have been the same if we omitted specifying the method argument since it’s the default method the function uses.
test_scores <- test_scores %>%
mutate(F1 = fscores(fit_2pl, method = "EAP"))
We can also calculate the model coefficient estimates using a generic function coef() which is used to extract model coefficients from objects returned by modeling functions. We will set the IRTpars argument to TRUE, which means slope intercept parameters will be converted into traditional IRT parameters.
coefs_2pl <- coef(fit_2pl, IRTpars = TRUE)
The resulting object coefs is a list, with one element for each question, and an additional GroupPars element that we won’t be using. The output is a bit long, so we’re only showing a few of the elements here:
coefs_2pl[1:3]
## $A1_B1
## a b g u
## par 1.1065245 -3.397220 0 1
## CI_2.5 0.9615469 -3.747013 NA NA
## CI_97.5 1.2515022 -3.047426 NA NA
##
## $A2_B0
## a b g u
## par 0.6021698 1.460565 0 1
## CI_2.5 0.5117190 1.237888 NA NA
## CI_97.5 0.6926207 1.683242 NA NA
##
## $A3_B2
## a b g u
## par 1.310811 -0.7101006 0 1
## CI_2.5 1.235152 -0.7592152 NA NA
## CI_97.5 1.386469 -0.6609860 NA NA
# coefs_2pl[35:37]
Let’s take a closer look at the first element:
coefs_2pl[1]
## $A1_B1
## a b g u
## par 1.1065245 -3.397220 0 1
## CI_2.5 0.9615469 -3.747013 NA NA
## CI_97.5 1.2515022 -3.047426 NA NA
In this output:
a is discriminationb is difficultyTo make this output a little more user friendly, we should tidy it such that we have a row per question. We’ll do this in two steps. First, write a function that tidies the output for one question, i.e. one list element. Then, map this function over the list of all questions, resulting in a data frame.
tidy_mirt_coefs <- function(x){
x %>%
# melt the list element
melt() %>%
# convert to a tibble
as_tibble() %>%
# convert factors to characters
mutate(across(where(is.factor), as.character)) %>%
# only focus on rows where X2 is a or b (discrimination or difficulty)
filter(X2 %in% c("a", "b")) %>%
# in X1, relabel par (parameter) as est (estimate)
mutate(X1 = if_else(X1 == "par", "est", X1)) %>%
# unite columns X2 and X1 into a new column called var separated by _
unite(X2, X1, col = "var", sep = "_") %>%
# turn into a wider data frame
pivot_wider(names_from = var, values_from = value)
}
Let’s see what this does to a single element in coefs:
tidy_mirt_coefs(coefs_2pl[1])
## # A tibble: 1 x 7
## L1 a_est a_CI_2.5 a_CI_97.5 b_est b_CI_2.5 b_CI_97.5
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 A1_B1 1.11 0.962 1.25 -3.40 -3.75 -3.05
And now let’s map it over all elements of coefs:
# use head(., -1) to remove the last element, `GroupPars`, which does not correspond to a question
tidy_2pl <- map_dfr(head(coefs_2pl, -1), tidy_mirt_coefs, .id = "Question")
A quick peek at the result:
tidy_2pl
## # A tibble: 39 x 7
## Question a_est a_CI_2.5 a_CI_97.5 b_est b_CI_2.5 b_CI_97.5
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 A1_B1 1.11 0.962 1.25 -3.40 -3.75 -3.05
## 2 A2_B0 0.602 0.512 0.693 1.46 1.24 1.68
## 3 A3_B2 1.31 1.24 1.39 -0.710 -0.759 -0.661
## 4 A4_B3 1.26 1.18 1.33 -0.586 -0.633 -0.538
## 5 A5_B4 1.04 0.960 1.11 -0.615 -0.684 -0.545
## 6 A6_B5 0.955 0.863 1.05 -2.28 -2.48 -2.09
## 7 A7_B6 1.43 1.34 1.51 -0.843 -0.892 -0.793
## 8 A8_B7 1.01 0.944 1.08 -0.541 -0.600 -0.482
## 9 A9_B8 1.63 1.53 1.74 -0.344 -0.390 -0.297
## 10 A10_B9 1.41 1.32 1.50 -0.691 -0.746 -0.635
## # ... with 29 more rows
And a nicely formatted table of the result:
gt(tidy_2pl) %>%
fmt_number(columns = contains("_"), decimals = 3) %>%
data_color(
columns = contains("a_"),
colors = scales::col_numeric(palette = c("Greens"), domain = NULL)
) %>%
data_color(
columns = contains("b_"),
colors = scales::col_numeric(palette = c("Blues"), domain = NULL)
) %>%
tab_spanner(label = "Discrimination", columns = contains("a_")) %>%
tab_spanner(label = "Difficulty", columns = contains("b_")) %>%
cols_label(
a_est = "Est.",
b_est = "Est.",
a_CI_2.5 = "2.5%",
b_CI_2.5 = "2.5%",
a_CI_97.5 = "97.5%",
b_CI_97.5 = "97.5%"
)
| Question | Discrimination | Difficulty | ||||
|---|---|---|---|---|---|---|
| Est. | 2.5% | 97.5% | Est. | 2.5% | 97.5% | |
| A1_B1 | 1.107 | 0.962 | 1.252 | −3.397 | −3.747 | −3.047 |
| A2_B0 | 0.602 | 0.512 | 0.693 | 1.461 | 1.238 | 1.683 |
| A3_B2 | 1.311 | 1.235 | 1.386 | −0.710 | −0.759 | −0.661 |
| A4_B3 | 1.255 | 1.183 | 1.327 | −0.586 | −0.633 | −0.538 |
| A5_B4 | 1.036 | 0.960 | 1.111 | −0.615 | −0.684 | −0.545 |
| A6_B5 | 0.955 | 0.863 | 1.047 | −2.283 | −2.475 | −2.090 |
| A7_B6 | 1.425 | 1.343 | 1.507 | −0.843 | −0.892 | −0.793 |
| A8_B7 | 1.010 | 0.944 | 1.076 | −0.541 | −0.600 | −0.482 |
| A9_B8 | 1.635 | 1.532 | 1.738 | −0.344 | −0.390 | −0.297 |
| A10_B9 | 1.406 | 1.316 | 1.496 | −0.691 | −0.746 | −0.635 |
| A11_B10 | 1.302 | 1.220 | 1.385 | −0.742 | −0.799 | −0.685 |
| A12_B0 | 0.783 | 0.673 | 0.893 | −1.705 | −1.936 | −1.473 |
| A13_B0 | 0.769 | 0.673 | 0.866 | 0.585 | 0.470 | 0.699 |
| A14_B12 | 1.361 | 1.281 | 1.440 | −0.500 | −0.547 | −0.453 |
| A15_B13 | 0.975 | 0.912 | 1.039 | −0.196 | −0.249 | −0.143 |
| A16_B14 | 1.249 | 1.167 | 1.332 | 1.200 | 1.129 | 1.271 |
| A17_B15 | 1.204 | 1.133 | 1.276 | −0.555 | −0.605 | −0.504 |
| A18_B16 | 0.879 | 0.815 | 0.944 | 0.601 | 0.534 | 0.668 |
| A19_B0 | 1.077 | 0.954 | 1.200 | 0.175 | 0.089 | 0.262 |
| A20_B17 | 1.807 | 1.699 | 1.914 | −0.957 | −1.004 | −0.909 |
| A21_B18 | 1.636 | 1.539 | 1.733 | −0.418 | −0.461 | −0.374 |
| A22_B19 | 1.636 | 1.536 | 1.735 | −0.515 | −0.561 | −0.469 |
| A23_B20 | 1.425 | 1.338 | 1.512 | −0.020 | −0.063 | 0.024 |
| A24_B21 | 1.227 | 1.149 | 1.305 | −0.715 | −0.773 | −0.658 |
| A25_B0 | 0.443 | 0.358 | 0.527 | −1.543 | −1.870 | −1.217 |
| A26_B22 | 1.631 | 1.486 | 1.776 | −2.000 | −2.124 | −1.877 |
| A27_B23 | 1.363 | 1.275 | 1.451 | −0.105 | −0.152 | −0.057 |
| A28_B24 | 1.292 | 1.205 | 1.380 | −0.339 | −0.393 | −0.284 |
| A29_B25 | 0.323 | 0.270 | 0.375 | −1.381 | −1.659 | −1.103 |
| A30_B26 | 0.956 | 0.890 | 1.023 | −0.899 | −0.973 | −0.826 |
| A31_B27 | 1.134 | 1.066 | 1.203 | 0.320 | 0.271 | 0.368 |
| A32_B28 | 1.175 | 1.094 | 1.257 | 1.191 | 1.119 | 1.262 |
| A33_B29 | 0.931 | 0.860 | 1.001 | −1.620 | −1.731 | −1.509 |
| A34_B30 | 0.548 | 0.482 | 0.614 | −2.447 | −2.738 | −2.156 |
| A35_B0 | 0.984 | 0.846 | 1.122 | −0.463 | −0.596 | −0.329 |
| A36_B0 | 0.928 | 0.802 | 1.053 | 0.662 | 0.542 | 0.783 |
| A0_B11 | 0.744 | 0.674 | 0.814 | −0.731 | −0.825 | −0.637 |
| A0_B31 | 0.748 | 0.670 | 0.825 | −0.068 | −0.155 | 0.020 |
| A0_B32 | 1.226 | 1.107 | 1.344 | 0.137 | 0.061 | 0.213 |
tidy_2pl %>%
mutate(qnum = parse_number(Question)) %>%
ggplot(aes(x = qnum, y = b_est, label = Question)) +
geom_errorbar(aes(ymin = b_CI_2.5, ymax = b_CI_97.5), width = 0.2) +
geom_text(colour = "grey") +
geom_point() +
theme_minimal() +
labs(x = "Question",
y = "Difficulty")
This shows the difficulty and discrimination of each of the questions, highlighting those that were added or removed:
tidy_2pl %>%
left_join(test_versions, by = c("Question" = "label")) %>%
mutate(qnum = parse_number(Question)) %>%
ggplot(aes(x = a_est, y = b_est, label = ifelse(outcome=="unchanged", "", Question), colour = outcome)) +
geom_errorbar(aes(ymin = b_CI_2.5, ymax = b_CI_97.5), width = 0.1, alpha = 0.5) +
geom_errorbar(aes(xmin = a_CI_2.5, xmax = a_CI_97.5), width = 0.1, alpha = 0.5) +
geom_text_repel() +
geom_point() +
theme_minimal() +
labs(x = "Discrimination",
y = "Difficulty")
Do students from different programmes of study have different distributions of ability?
Compare the distribution of abilities in the year groups (though in this case there is only one).
ggplot(test_scores, aes(F1, y = year, fill = as.factor(year), colour = as.factor(year))) +
geom_density_ridges(alpha=0.5) +
scale_x_continuous(limits = c(-3.5,3.5)) +
labs(title = "Density plot",
subtitle = "Ability grouped by year of taking the test",
x = "Ability", y = "Density",
fill = "Year", colour = "Year") +
theme_minimal()
## Picking joint bandwidth of 0.182
plot(fit_2pl, type = "infoSE", main = "Test information")
plot(fit_2pl, type = "infotrace", main = "Item information curves")
plot(fit_2pl, type = "score", auto.key = FALSE)
We can get individual item surface and information plots using the itemplot() function from the mirt package, e.g.
mirt::itemplot(fit_2pl, item = 1,
main = "Trace lines for item 1")
We can also get the plots for all trace lines, one facet per plot.
plot(fit_2pl, type = "trace", auto.key = FALSE)
Or all of them overlaid in one plot.
plot(fit_2pl, type = "trace", facet_items=FALSE)
An alternative approach is using ggplot2 and plotly to add interactivity to make it easier to identify items.
# store the object
plt <- plot(fit_2pl, type = "trace", facet_items = FALSE)
# the data we need is in panel.args
# TODO - I had to change the [[1]] to [[2]] since the plt has two panels for some reason, with the one we want being the 2nd panel
plt_data <- tibble(
x = plt$panel.args[[2]]$x,
y = plt$panel.args[[2]]$y,
subscripts = plt$panel.args[[2]]$subscripts,
item = rep(colnames(item_scores), each = 200)
) %>%
mutate(
item_no = str_remove(item, "A") %>% as.numeric(),
item = fct_reorder(item, item_no)
)
head(plt_data)
## # A tibble: 6 x 5
## x y subscripts item item_no
## <dbl> <dbl> <int> <fct> <dbl>
## 1 -6 0.0531 201 A1_B1 NA
## 2 -5.94 0.0566 202 A1_B1 NA
## 3 -5.88 0.0603 203 A1_B1 NA
## 4 -5.82 0.0642 204 A1_B1 NA
## 5 -5.76 0.0683 205 A1_B1 NA
## 6 -5.70 0.0727 206 A1_B1 NA
plt_gg <- ggplot(plt_data, aes(x, y,
colour = item,
text = item)) +
geom_line() +
labs(
title = "2PL - Trace lines",
#x = expression(theta),
x = "theta",
#y = expression(P(theta)),
y = "P(theta)",
colour = "Item"
) +
theme_minimal() +
theme(legend.position = "none")
ggplotly(plt_gg, tooltip = "text")
knitr::knit_exit()